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O ABSTRACT 

I ■ We present a new evolutionary model that describes the population properties of ra- 

dio sources at frequencies ;$ 5 GHz, thus complementing the De Zotti et al. (2005) 

model, holding at higher frequencies. We find that simple analytic luminosity evolu- 
tion is still sufficient to fit the wealth of available data on local luminosity functions, 
multi-frequency source counts, and redshift distributions. However, the fit requires 
a luminosity-dependent decline of source luminosities at high redshifts, at least for 
steep-spectrum sources, thus confirming earlier indications of a "downsizing" also for 
radio sources. The upturn of source counts at sub-mJy levels is accounted for by a 
\^ | straightforward extrapolation, using the empirical far-IR/radio correlation, of evolu- 

, tionary models matching the far-IR counts and redshift distributions of star-forming 

galaxies. We also discuss the implications of the new model for the interpretation 
of data on large-scale clustering of radio sources and on the Integrated Sachs- Wolfe 
(ISW) effect, and for the investigation of the contribution of discrete sources to the 
' extragalactic background. As for the ISW effect, a new analysis exploiting a very clean 

CMB map, yields at a substantially higher significance than reported before. 
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1 INTRODUCTION (Best, Rottgering & Lehnert 1999; Willott et al. 2002; Jack- 

. , . ,. son et al. 2002; Brookes et al. 2008; Gendre & Wall 2008). 

1 he still most widely used evolutionary models lor radio „ , . , , . a . r , . , , , . „ 

. — ^ „„„„, , _ Complete redshift information on source samples at drfter- 

sources at c ^ 5 GHz date back to the 1990 s Dunlop & , n ., . . ., . , , , , ,. 

v ent faux limits is especially critical to clarity the evolutionary 

Peacock 1990; Toffolatti et al. 1998). Although these mod- properties. For example, as shown by Brookes et al. (2008), 

els proved to be very useful benchmarks even today, several the observed redshift distribution of the CENSORS sample, 

of their building blocks need to be updated, starting from complete down to S1.4G&, = 7.2 mjy, is not well reproduced 

the underlying cosmology (a matter dominated flat universe by any of the Dunlop & Peacock's (1990) models. All that 

with H = 50 kms -1 Mpc -1 ). Even more importantly, large points to the need of updated models. 

amounts of new data have been accumulating since then, Many kgy sdentific issues in this field are stm open: 

including new deep/large area surveys (see De Zotti et al. what ig the nature q{ evolution? is it i uminosity dep en- 

2009 for a review and references) , identifications of sub-mJy dent? ig there a ^ evidence for a decline of the comoving 

sources at 1.4 GHz (Seymour et al. 2008; Bondi et al. 2008; space dengity of radio sources ftt high redsMfts? are flat _ and 

Ibar et al. 2009; and references therein), accurate deter- ste ep-spectrum sources obeying the same evolution laws, as 

minations of the local luminosity function (Magliocchetti expected in the framework of 'unification' schemes? which 

et al. 2002; Condon et al. 2002; Best et al. 2005; Mauch sourceg afe responsible for the upturn of differential source 

& Sadler 2007), redshift distributions for complete samples counts at sub-mJy levels? 

In this paper we will investigate to what extent the new 
data allow us to progress towards answering these questions. 

* E-mail: marcella.massardi@oapd.inaf.it The analysis is similar and complementary to the one by De 
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Zotti et al. (2005), holding for v ^ 5 GHz. In fact, man- 
ageable models cannot fully take into account the complex 
and multiform spectral behaviour of radio sources. They 
must adopt schematic descriptions, that most frequently 
boil down to considering two populations, steep- and flat- 
spectrum sources, each with a single power-law spectrum, 
S oc v~ a . At best, a Gaussian distribution of spectral in- 
dices is allowed for. In contrast, real source spectra, espe- 
cially those classified as "flat", are generally not power-laws, 
but have complex shapes, showing spectral bumps, Batten- 
ings or inversions (i.e. flux increasing with increasing fre- 
quency) , frequently bending to steeper power-laws at higher 
frequencies. This means that power-law approximations can 
only hold for limited frequency ranges. In particular, the De 
Zotti et al. (2005) model becomes increasingly inaccurate 
with decreasing frequency below 5 GHz. Furthermore, the 
high frequency data sets analyzed by De Zotti et al. (2005) 
are dominated by compact, flat-spectrum sources and do 
not strongly constrain the evolutionary properties of the 
steep-spectrum population. Much stronger constraints are 
provided by the very abundant data collected at low fre- 
quencies, which are the subject of the present analysis. 

In §[2] we present our parametrization of the evolving 
luminosity functions of AGN-powered radio sources and its 
application to estimate source counts and redshift distribu- 
tions. In §[3] we deal with the contributions of star-forming 
galaxies to the source counts. In §[3] the integrated emis- 
sion of sources is compared to observational estimates of 
the radio background. In §[5] we discuss the implications of 
the new model for the interpretation of the angular correla- 
tion function of radio sources and of data on the Integrated 
Sachs-Wolfe effect. Finally, in §[6] we summarize our conclu- 
sions. 

We have adopted a flat ACDM cosmology with A = 0.7 
and H = 70kms _1 Mpc -1 . 



2 AGN POWERED SOURCES 
2.1 Data sets 

2.1.1 Local luminosity functions 

A spectacular progress towards the determination of the lo- 
cal radio luminosity function (LRLF) at 1.4 GHz has been 
achieved combining large-area spectroscopic surveys (Las 
Campanas, SDSS, 2dF and 6dF) with the NVSS and FIRST 
surveys (Machalski & Godlowski 2000; Magliocchetti et al. 
2002; Sadler et al. 2002; Best et al. 2005; Mauch & Sadler 
2007). Condon et al. (2002) estimated the LRLF by cross- 
matching NVSS sources brighter than 2.5 mjy with galaxies 
brighter than m p = 14.5 mag in the Uppsala Galaxy Cata- 
logue (UGC; Nilson 1973). 

We will mostly refer to the analysis by Mauch & Sadler 
(2007) who derived the LRLF from a large, well-defined sub- 
sample of radio sources with Si.4GHz ^ 2.8 mjy from the 
NVSS catalogue, identified with galaxies observed in the 
Second Incremental Data Release (Jones et al. 2005) of the 6 
degree Field Galaxy Survey (6dFGS). It may be noted, how- 
ever, that the NVSS is only 50% complete at 2.5 mjy, and its 
completeness rises rapidly to 99% at 3.4 mjy (Condon et al. 
1998). The incompleteness is therefore still quite significant 
(~ 30%) at 2.8 mjy. This likely implies that the true uncer- 



tainties are larger than the quoted statistical errors. In fact, 
estimates of the local luminosity functions using the samples 
with complete (or almost complete) redshift information (see 
§ !2.1.3p . and in particular the steep-spectrum sub-sample of 
Kiihr et al. (1981) and the CoNFIG sample (Gendre & Wall 
(2008) give, on average, local space densities about 0.2 dex 
lower in the luminosity range for which these samples have 
sufficient statistics (around log L\a GHz(ergs _1 Hz -1 ) ~ 32; 
see also Rigby et al. 2008). In our analysis we have therefore 
adopted an error of at least 0.2 dex for the Mauch & Sadler 
(2007) LRLF of radio AGNs. Also, Mauch & Sadler (2007) 
have included in their "local" sample sources up to z ~ 0.2. 
However, the V/V ma x test for the AGN sub-sample showed 
evidence for cosmological evolution at more than 3cr. We 
have therefore restricted ourselves to sources at z ^ 0.05, i.e. 
to AGNs with log L1.4 GHz(erg s _1 Hz _1 ) < 33.3, for which 
evolution is expected to have a negligible effect. 

A further important recent advance is the derivation 
of separate luminosity functions for the star-forming galax- 
ies and for the AGNs (see the upper panel of Fig. [TJ). This 
has allowed us to easily deal with the different evolutionary 
properties of the two source populations. Note that star- 
forming galaxies have, in most cases, steep spectra. 

On the other hand, recent studies did not attempt to 
disentangle the contributions to the LRLF of sources of dif- 
ferent spectral classes. The time-honoured determinations 
of the LRLFs of flat- and steep-spectrum sources by Dun- 
lop & Peacock (1990) and Toffolatti et al. (1987) are com- 
pared, in the top panel of Fig. [T] with the total LRLF by 
Mauch & Sadler (2007). The agreement with the Toffolatti 
et al. (1987) estimates is extremely good, while the estimate 
by Dunlop & Peacock (1990) is somewhat high at low lu- 
minosities. As shown in the central panel of the figure, for 
log L[erg s _1 Hz _1 ] £ 30, the global LRLF of steep-spectrum 

sources closely matches Mauch & Sadler's (2007) determi- 
nation of the starburst LRLF, while at higher luminosities 
it matches the AGN LRLF. The bottom panel shows that 
the Toffolatti et al. (1987) LRLF of flat-spectrum sources 
joins smoothly with the estimated LRLF of BL Lac objects 
obtained by Padovani et al. (2007). We have included the 
Toffolatti et al. (1987) LRLFs of flat- and steep-spectrum 
sources and the Padovani et al. (2007) LRLF for BL Lacs in 
our data-base. 

Donoso et al. (2009) presented a determination of the 
1.4 GHz radio luminosity function at z ~ 0.55 by cross corre- 
lating the NVSS and FIRST catalogs with the MegaZ Lumi- 
nous Red Galaxy (MegaZ-LRG) catalog derived from Sloan 
Digital Sky Survey imaging data. The combined catalog has 
an estimated reliability of ~ 98.3 per cent and completeness 
level of about 95 per cent. Accurate photometric redshifts 
are available for all the objects. Although the sample is very 
large, and the statistical errors on the luminosity function 
correspondingly small, it is important to note that the radio 
luminosity function misses radio galaxies falling outside the 
MegaZ-LRG colour selection criteria, as well as radio-loud 
quasars. On the whole, the Donoso et al. estimate may be 
too low by 0.1-0.2 dex (see Fig. [2}. For this reason we have 
adopted a positive error of at least 0.2 dex for this dataset. 
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Figure 1. Comparison of model with observed local radio luminosity functions (LRLFs) at 1.4 GHz. The upper panel refers to the global 
(flat- plus steep-spectrum source, FS and SS respectively) luminosity function. The solid line shows the best fit LRLF of radio AGNs 
and the dotted and the dashed lines give the contributions of flat- and steep-spectrum components, respectively. The (triple) dot-dashed 
line shows the LRLF of star-forming (SF) galaxies, obtained by (linearly) non-linearly extrapolating the 60 jim luminosity function by 
means of the well established far-IR/radio correlation (see Q. The estimate by Toffolatti et al. (1987; triangles) is in remarkably good 
agreement with the the recent assessment by Mauch & Sadler (2007; diamonds), over the whole luminosity range. The estimate by 
Dunlop & Peacock (1990, their "determination 1"; squares) agrees well with for logL(ergs _1 Hz ) > 29.5, but is systematically higher 

at lower luminosities. Error bars are omitted to avoid overcrowding. The central panel refers to steep-spectrum sources. The diamonds 
and the asterisks show, respectively, the LRLFs of radio AGNs and of star-forming galaxies derived by Mauch & Sadler (2007) who did 
not distinguish between steep- and flat-spectrum sources. The triangles and the squares show the estimates by Toffolatti et al. (1987) 
and Dunlop & Peacock (1990, "determination 1") for steep spectrum sources only, but including both AGNs and star-forming galaxies. 
The solid line is the best fit LF of steep-spectrum AGNs and the dot-dashed line that of star- forming galaxies. The latter is the result of 
the non- linear extrapolation of the local 60 fim luminosity function (see text). Finally, the bottom panel shows the LRLF of flat-spectrum 
sources (AGNs only). The observational estimates are, again, from Toffolatti et al. (1987) and Dunlop &; Peacock (1990). At the highest 
luminosities these join smoothly with the LRLF of BL Lacs by Padovani et al. (2007; asterisks). At the lowest luminosities the LRLF of 
Dunlop & Peacock (1990, "determination 1"; squares) is slightly high in comparison with both Toffolatti et al. (triangles) and the global 
(flat- plus steep-spectrum) LRLF of AGNs by Mauch & Sadler (2007; diamonds), shown for comparison. The dashed and the dotted 
lines are the contributions of BL Lacs and of Flat Spectrum Radio Quasars (FSRQs), respectively, to the best fit LRLF of flat-spectrum 
sources (solid line). The Dunlop & Peacock (1990) and Toffolatti et al. (1987) data have been rescaled to the value of Ho used in this 
paper and extrapolated from 2.7 to 1.4 GHz using a spectral index of 0.8 for steep-spectrum sources, and of 0.1 for flat-spectrum sources. 



2.1.2 Source counts 

An updated compilation of radio source counts is presented 
by De Zotti et al. 

(20090- 

Our analysis deals with those in 
the frequency range 0.15 ^ v ^ 5 GHz. While above ~ 1 mjy 



1 Tables of the source counts and model re- 
sults can be downloaded from the webpage 
http: / /web. oapd.inaf.it / rstools /sreent /srccnt.tables.html 



counts are dominated by AGNs, the contribution of star- 
forming galaxies becomes increasingly important at fainter 
flux levels (Windhorst et al. 1985; Benn et al. 1993; Hopkins 
et al. 1998; Seymour et al. 2004; Muxlow et al. 2005; Moss et 
al. 2007). Therefore sub-mjy source counts were not taken 
into account in our model fitting process, unless the contri- 
butions of radio AGNs were reliably singled out. This is the 
case for Seymour et al. (2008) who used several diagnostics 
to determine the source counts of AGNs down to less than 
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Figure 2. Model 1.4 GHz RLF at different redshifts, specified in the inset, for FSRQs (upper panel) and steep-spectrum (SS) AGNs 
(lower panel). The Donoso et al. (2009) 1.4 GHz RLF at z ~ 0.55 is also shown (asterisks) for comparison (the model RLF at the 
same redshift is given by the solid line). This data does not distinguish between flat and steep-spectrum objects, but at 1.4 GHz and 
at z ~ 0.55 the sample is expected to be mostly constituted by steep-spectrum sources; it is therefore compared with the model for SS 
AGNs. 



0.1 mJy at 1.4 GHz; their faint AGN counts were included in 
the dataset used for the fitting. As shown in §[3j the contri- 
bution to the counts of star-forming galaxies can be obtained 
by straightforwardly extrapolating models fitting mid- and 
far-IR data, without introducing any adjustable parameter. 
On the bright side, an issue of concern is the possibility that 
giant radio sources, which have global flux density above the 
selection limit of a sample are missed because of resolution 
effects (Riley 1989). The observed counts may therefore be 
lower limits to the true counts. 

In the case of small-area surveys, the sampling variance 
due to source clustering can be a further significant addi- 
tion to the uncertainties. The total variance of the counts is 
(Peebles 1980); 



(n) 



(n) 



i 2 



with 



w{9) dSli dtt 2 



(1) 



(2) 



where 9 is the angle between the solid angle elements dQ\ 
and dQ2, and the integrals are over the solid angle covered 



by the survey. The angular correlation function of NVSS and 
FIRST sources is consistent with a power-law shape (Blake 
& Wall 2002; Overzier et al. 2003): 



to (61) ~ lO" 3 (0/deg)- 



(3) 



for angular separations up to at least 4° . Inserting eq. ([3]) in 
eq. we get 



a 2 = 2.36 x 10~ 3 (n/deg 2 )" 



(4) 



This contribution was summed in quadrature to the Poisson 
errors for surveys over areas ^ 25deg 2 . However the actual 
differences between published source counts from indepen- 
dent fields are much larger than those accounted for by Pois- 
son plus clustering fluctuations (see Fig. 4 of Condon 2007). 
As discussed by Condon (2007) the dominant contributions 
to such differences are likely due to calibration errors, resolu- 
tion corrections, corrections for primary-beam attenuation, 
applied in different ways by different observers. For example, 
the deepest surveys cover a single primary-beam area with 
nonuniform sensitivity, so the counts contain large correc- 
tions for primary gain; and surveys made with high angular 
resolution must be corrected for missing extended sources. 
This means that using the Poisson plus clustering errors to 
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Figure 3. Comparison of our best fit model with the observed source counts (multiplied by a factor c for a better visualization) at 
several frequencies. From top to bottom: 150 MHz (c = 10000), 325 MHz (c = 1000), 408 MHz (c = 100), 610 MHz (c = 10), 1.4 GHz 
(c = 1), 5.0 GHz (c = 0.1). The solid and the dotted lines show our model counts for AGNs and starforming galaxies, respectively, while 
the dashed lines show the sum of the two contributions. The dot-dashed line at 5 GHz shows the counts yielded by the De Zotti et al. 
(2005) model for AGNs; they are essentially indistinguishable from those of the present model, except at sub-mjy levels. For a complete 
list of references and tabulations of the data points see De Zotti et al. (2009). 



weight the source count data in the comparison with the 
model may be misleading. On the other hand, it is diffi- 
cult to find consistent alternative weighting schemes, and 
this leaves some level of discretion in the comparison of the 
model to the data. 

2.1.3 Redshift distributions 

As mentioned in §[T] particularly important new constraints 
on evolutionary models have been provided by the redshift 
distribution for the Combined EIS-NVSS Survey Of Radio 
Sources (CENSORS; Brookes et al. 2008). We have consid- 
ered the complete sample of 136 AGNs with Sii m ,i.4GHz = 
7.2 mjy over a 6 deg 2 area. The 9 sources having a lower 
limit for the redshift have been taken into account using 
the Kaplan-Meier estimator [routine KMESTM in the soft- 
ware package ASURV Rev 1.2 (Isobe & Feigelson, 1990), 
which implements the survival analysis methods presented 
in Feigelson & Nelson (1985) and Isobe, Feigelson & Nelson 
(1986)]. Poisson errors have been adopted. 

While the CENSORS sample is unique in providing in- 
formation on the redshift distribution of faint sources, spec- 



troscopically complete (or close to completeness) samples of 
bright sources include: 

• the 3CRR sample (Laing, Riley & Longair 1983) con- 
taining a total of 173 sources, complete for Si78MHz > 
10.9 Jy and with spectral index information and redshifts 
for all sources 

• the Best, Rottgering, & Lehnert (1999) sample of 178 
sources brighter than AosmHz = 5 Jy, with complete redshift 
information (Best et al. 2003) 

• the Wall & Peacock (1985) 2-Jy sample of 239 sources 
selected at 2.7 GHz, only 5 of which are missing a redshift 
estimate 

• the 7C redshift survey (7CRS; Willott et al. 2002; Lacy, 
Bunker & Ridgway 2000): 103 sources with Sisimhz ^ 0.5 Jy 
and full redshift information. 

• the Combined NVSS-FIRST Galaxies (CoNFIG) sam- 
ple (Gendre & Wall 2008) comprising 274 sources with 
5i.4GHz > 1.3 Jy. Spectroscopic redshifts are available for 
230 sources (84%) and photometric estimates were made for 
14 sources, bringing the total redshift coverage to 89%. 

While all these samples do not distinguish between steep- 
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Figure 4. Contributions of different source populations (see inset) to the source counts at 0.5, 1.4, and 5.0 GHz. 



and fiat-spectrum sources, the Kiihr et al. (1981) catalog 
provides also spectral information. The sample comprises 
518 sources to a 5 GHz flux density limit of 1 Jy, over 
an area of 9.811 sr. Based on the catalogued spectral in- 
dices, 299 sources are flat-spectrum (a < 0.5, S oc v~ a ) 
and 219 are steep-spectrum (a ^ 0.5). The former pop- 
ulation includes 212 flat-spectrum radio quasars (FSRQs), 
200 of which (94%) have measured redshift, 26 BL Lacs (20 
of which, 77%, with measured redshift, and 61 either clas- 
sified as galaxies or missing a morphological classification. 
We have distributed the latter sources among FSRQs and 
BL Lacs in proportion of the fraction of classified sources. 
The fraction of steep-spectrum sources with measured red- 
shift is of 81%. 

Jackson et al. (2002) presented 99% complete optical 
identifications of the 878 sources of the Parkes 2.7 GHz 
Quarter- Jansky flat-spectrum sample. Redshifts are avail- 
able for 58% of the sources. Following De Zotti et al. (2005) 
we have denned a complete sub-sample aiming at maximiz- 
ing the fraction of sources with measured redshift. The sub- 
sample includes the sources brighter than 0.25 Jy detected 
in areas surveyed at least down to this limit at 5 > —50°, 
for a total of 514 objects of which 370 are FSRQs (93% with 
redshift), 47 are BL Lacs (21% with redshift), 95 are fiat- 
spectrum radio galaxies, and 2 are unidentified. Only for 
FSRQs a meaningful redshift distribution can be derived. 



2.2 The model 

Our model includes two flat-spectrum populations with dif- 
ferent evolutionary properties (FSRQs, and BL Lacs) and 
a single steep-spectrum population. For each population 
we adopt a simple power-law spectrum, S oc v~ a , with 

QFSRQ = QBLLac = 0.1, and Qfstcop = 0.8. 

The epoch-dependent comoving luminosity functions (in 
units of Mpc -3 (dlogL) -1 ) are modeled as double power- 
laws: 



<£>(£, z) = 



n 



(5) 



The high-z evolution of radio sources has been the subject 
of a hot, still not completely settled, debate (Dunlop & Pea- 
cock 1990; Shaver et al. 1996; Jarvis & Rawlings 2000; Wall 
et al. 2005; Cruz et al. 2007). Standard luminosity evolu- 
tion models (£*(z) = L*(0) exp[k evo r(z)], r(z) being the 
lookback time in units of the Hubble time, Hq 1 ), without 
a high-z decline of the comoving luminosity function, did 
not yield any acceptable fit to the data on steep-spectrum 
AGNs and on FSRQs. For these source populations we have 
therefore adopted a more complex luminosity evolution law: 



L*(z) = L*(0) dex ■ 



Stop Z 



(l-m ev ) 

nop 



'/(l + mev))],(6) 



where z top is the redshift at which L*(z)/L*(0) reaches its 
maximum and |m ov | < 1. 
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Figure 5. Normalized rcdshift distributions yielded by our best fitting model (solid lines) compared with data by Gendre & Wall (2008), 
3CRR (Laing, Riley & Longair 1983), Brookes et al. (2008), Best, Rottgering, & Lehnert (1999), Wall & Peacock (1985), 7CRS (Willott 
ct al. 2002; Lacy, Bunker & Ridgway 2000), Jackson et al. (2002; quarter-Jy, QJy, sample, which includes only FSRQs) and the rcdshift 
distributions for the steep-spectrum, FSRQ, and BL Lac sub-samples of the Kuhr et al. (1981) sample; see insets. 
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Figure 6. Contributions of the different populations to the redshift distributions yielded by our best fitting model at different frequencies 
(0.5, 1.4, 5.0 GHz, from top to the bottom) and different flux limits (0.001, 0.01, 1 Jy from left to right). 
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Table 1. Best fitting model parameters. The spectral indices were 



fixed at Qfsrq = QBLLac = 0.1, and a ste 



: Q.f 



Table 2. Contributions to \ 2 of each data set for the best fit 
values of the parameters (see Table [TJ. 



Parameters 


FSRQ 


BLLac 


SS-AGNs 


a 


0.760 


0.723 


0.559 


b 


2.508 


1.618 


2.261 


log n 


-10.382 


-6.879 


-5.970 


log MO) 


34.323 


32.638 


32.490 


^evo 


-0.996 


0.208 


1.226 


Ztop.O 


1.882 


1.282 


0.977 




0.018 




0.842 




-0.166 


1 


0.262 



Although this improved substantially the quality of the 
fit, significant discrepancies between the model and the data 
remained. In particular, the model fell short at reproduc- 
ing the highest z points of the redshift distributions of the 
Brookes et al. (2008) sample and the Kiihr et al. (1981) 
steep-spectrum sub-sample (see Fig. [5]). We can force the 
model to fit those points, but at the cost of visibly wors- 
ening the fit to the counts and/or to the local luminosity 
functions. A better global fit was achieved allowing z top to 
vary with luminosity as: 



Ztop = Ztop,0 + 



1 + L*(0)/L' 



(7) 



ztop.o and Sztop being free parameters. The fit gave posi- 
tive values for both of them, consistent with earlier indica- 
tions that the high-z decline of the space density is more 
pronounced and starts at lower redshifts for less powerful 
sources (Hook et al. 1998; Waddington et al. 2001; Vigotti 
et al. 2003; Cirasuolo et al. 2005, 2006; Wall et al. 2009 in 
preparation). The luminosity dependence of the high-z de- 
cline is qualitatively similar to the downsizing observed for 
galaxies and optically and X-ray selected quasars (Cowie et 
al. 1996; Barger et al. 2005; Perez-Gonzalez et al. 2008). 

Actually, as shown by Table [T] the luminosity depen- 
dence of the peak redshift required by the data is substantial 
for the steep-spectrum population but is weak (8z t op "C 1) 
for FSRQs. This is expected since the contributions to the 
observables (source counts and redshift distributions) of 
moderate to low luminosity flat-spectrum sources are dom- 
inated by BL Lacs, so that the evolution of low luminosity 
FSRQs is poorly constrained. 

In the case of BL Lacs we have adopted an evolutionary 
form similar to those used for the FSRQ and steep-spectrum 
sources, but, since the data are not enough to constrain sev- 
eral evolutionary parameters for this population, we have 
used the simpler form of eq. (JUJ), with m cv = 1 (fixed) and 
a luminosity-independent redshift cut-off. 



2.3 Multidimensional minimization algorithms 

The best-fit model has been found by minimizing the sum 
of x 2 values obtained comparing the model with data 
on luminosity functions, redshift distributions and source 
counts. The minimization was first performed using the 
multi-dimensional downhill simplex method (exploiting the 
AMOEBA routine, Press et al. 1992). The method is compu- 
tationally fast even in the case of a large number of param- 
eters. However, it does not allow us to bracket a solution, 



Data set 


x 2 


Ndata 


X 2 / N dat a 


150 MHz counts [1] 


110.0 


45 


2.444 


325 MHz counts [1] 


52.2 


17 


3.071 


408 MHz counts [1] 


57.7 


41 


1.311 


610 MHz counts [1] 


63.8 


37 


1.724 


1.4 GHz counts [1] 


363.0 


115 


3.156 


4.8 GHz counts [1] 


50.7 


16 


1.102 


z distr 121 


29.1 


20 


1.455 


z distr. [3] FSRQ 


24.9 


16 


1.556 


z distr. [4] FSRQ 


9.5 


21 


0.396 


z distr. [4] BLLac 


12.5 


19 


0.658 


z distr. [4] SS 


36.8 


35 


1.051 


z distr. [5] 


23.0 


26 


0.885 


z distr. [6] 


11.8 


21 


0.561 


z distr. [7] 


26.7 


11 


0.651 


z distr. [8] 


32.6 


26 


1.254 


z distr. [9] 


22.2 


29 


0.766 


LLF [10] AGN 


7.4 


12 


0.617 


LLF [11] SS 


13.6 


5 


2.720 


LLF [12] BLLac 


33.7 


12 


2.808 


LF [13] at z = 0.55 


1.4 


8 


0.553 


Total 


963.0 


598 


1.610 



References: [1] De Zotti et al. (2009); [2] Brookes et al. (2008); 
[3] Jackson et al. (2002); [4] Kiihr et al. (1981); [5] Gendre & 
Wall (2008); [6] Laing, Riley & Longair (1983); [7] Willott et al. 
(2002); [8] Best, Rottgering, & Lehnert (1999); [9] Wall & Peacock 
(1985); [10] Mauch & Sadler (2007); [11] Toffolatti et al. (1987); 
[12] Padovani et al. (2007); [13] Donoso et al. (2009) 



but, starting from an initial guess, modifies the parameters 
defining its path towards the closest minimum. So, if the 
parameter space has several local minima, the method can 
be trapped at a secondary minimum or can find a minimum 
for physically implausible values of some parameters. 

In order to check the sensitivity of the downhill simplex 
method to the initial guess and to confirm the final results, 
we have also exploited the simulated annealing minimiza- 
tion technique (Aarts & Korst 1989). This method exploits 
an analogy between the way in which a metal cools and 
freezes into a minimum energy crystalline structure and the 
search for a minimum in a more general system. The ma- 
jor advantage of simulated annealing over other methods is 
its ability to avoid becoming trapped at secondary minima. 
Also it allows us to specify the plausible range for each pa- 
rameter. The algorithm employs a random search which not 
only accepts changes that decrease the merit function, but 
also some changes that increase it. The probability of the lat- 
ter to occur is inversely proportional to the increase in the 
merit function. The disadvantage is that it is computation- 
ally demanding, especially for a large number of parameters, 
as in our case. Therefore, we have used it only to check the 
stability of the results obtained with the other method. 

Care must be taken to avoid that some data sets influ- 
ence too heavily the results. In particular the redshift dis- 
tributions, which contain far fewer points than the counts, 
and have larger error bars may be under-weighted. To check 
whether this may seriously affect the results we have re- 
peated the minimization increasing the weights of the red- 
shift distributions by up to a factor of 10. Since the best 
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fit values of the parameters turned out to be remarkably 
stable we have adopted those obtained without introducing 
additional weights. 

The best fit values of the parameters are listed in Ta- 
ble [T] and Table [2] details the contributions to \ 2 °f each 
data set. The fact that fits of the source counts yield values 
of xV^data as high as 2 or even 3 is not really worrisome 
since large contributions to \ 2 come from data points from 
different surveys which are inconsistent with other data; as 
noted in § 12.1,21 differences among source counts from inde- 
pendent fields are much larger than the statistical errors. 
Allowing for that, the quality of the fit appears to be ac- 
ceptable for all data sets. 

Comparisons of the best fit model with several data sets 
are shown in Figs. \T\ [3] and The disaggregation of the 
counts yielded by the best fit model into the contributions 
of the different populations is shown in Fig. [4] The contribu- 
tion of different populations to the redshift distributions for 
different flux density limits at different frequencies is shown 
in Fig. [6] The plotted redshift distributions are normalized, 
i.e. the number of sources in each redshift bin is divided by 
the total number of sources in the sample. The contribu- 
tion of star-forming galaxies to the sub-mjy counts will be 
discussed in the next section. 

As mentioned in § 12.21 we have tried different evolu- 
tionary models. Simple standard luminosity evolution mod- 
els (L*(z) = L*(0) exp[fc cvo r(z)], t(z) being the lookback 
time in units of the Hubble time) produced absolutely un- 
acceptable fits for both steep- and flat-spectrum sources. 
The fit strongly improved when the evolution function of 
both steep-spectrum sources and FSRQs included a decline 
at high-z, described by eq. (j6j), while keeping simple lumi- 
nosity evolution for BL Lacs. The global minimum \ 2 be- 
came 1140 (to be compared with Xmin = 963 obtained for 
the model in Table Q] see Table [2}. A further improvement 
(Xmin — 1060) was achieved allowing z top to decrease with 
decreasing luminosity, according to eq. (J7J). Actually, the de- 
crease of Xmin was entirely due to data on steep-spectrum 
sources; setting Sz t0 p = for flat-spectrum sources did not 
change in any appreciable way the quality of the fit. Finally, 
we used also for BL Lacs the evolutionary law of eq. ([6]); this 
also improved the fit, yielding our final Xmin — 963. 



3 STAR-FORMING GALAXIES 

The radio emission of star-forming galaxies correlates with 
their star formation rate, as demonstrated by the well- 
established tight correlation with far-IR emission (Helou et 
al. 1985; Gavazzi et al. 1986; Condon 1992; Garrett 2002). 
Yun et al. (2001) found that the overall trend in the range 
L„(60^m) ~ 10 30 -10 32 - 5 ergs -1 Hz -1 is indistinguishable 
from a linear relation: 

Liagbx = 1.16 x 10~ 2 L„(60^m). (8) 

Galaxies with L„(60/im) < 10 30 ergs -1 Hz -1 are found to 
have radio to far-IR luminosity ratios systematically lower 
than those given by eq. ([8]). The apparent deviation from 
linearity in the radio/far-IR correlation at low luminosities 
is supported by a comparison of 60 /im and 1.4 GHz local lu- 
minosity functions (Yun et al. 2001; Best et al. 2005). Simply 
shifting the 60 fim luminosity function (Saunders et al. 1990; 



Takeuchi et al. 2003) along the luminosity axis according to 
eq. iJSJ yields a good match to the radio luminosity func- 
tion (Best et al. 2005; Mauch & Sadler 2007, see Fig. for 
Liaghz Si 10 28 ' 5 ergs -1 Hz -1 . At lower luminosities, how- 
ever, the extrapolated luminosity function lies increasingly 
above the observed one (Fig. [TJl . Full agreement is recovered 
by replacing eq. © with 

1.16 x 10- 2 L b 

Li 4GHz q— i 1 (yj 

(L b /L,(60fim)) 31 + L h /L v (<BOfxm) 

in which L b = 8.8 x 10 29 ergs" 1 Hz" 1 . 

The tight empirical relationship between radio and far- 
IR luminosities for star-forming galaxies allows us to take 
advantage of the wealth of data at far-IR/sub-mm wave- 
lengths to derive the radio evolution properties. We expect a 
different evolution for starburst and normal late-type galax- 
ies as the starburst activity is likely triggered by interactions 
and mergers that were more frequent in the past, while in 
normal galaxies the star-formation rate has probably not 
changed much over their lifetimes. The bulk of the sub-mm 
counts measured by SCUBA surveys (Scott et al. 2006; Cop- 
pin et al. 2006) is due to yet another population, the sub-mm 
galaxies (SMGs), interpreted as proto-spheroidal galaxies in 
the process of forming most of their stars (Granato et al. 
2004). 

A straightforward extrapolation to radio frequencies of 
the evolutionary models by Negrello et al. (2007) for the 
three populations (normal, starburst and sub-mm galax- 
ies), exploiting eq. © and the SEDs of NGC6946 for nor- 
mal late-type galaxies and of Arp220 for starburst and 
proto-spheroidal galaxies, yields the counts shown in Fig. [3] 
Adding the contributions of star-forming galaxies to the 
counts yielded by the model for radio AGNs described in 
ij |2.2l the sub-mjy counts are nicely reproduced. The counts 
at tens of /iJy levels turn out to be dominated by star- 
forming galaxies. Figure [S] details the contribution of each 
population of star-forming galaxies to the (sub)mjy counts 
at 1.4 GHz; the counts of radio AGNs are also shown for 
comparison. 



4 CONTRIBUTIONS TO THE RADIO 
BACKGROUND 

The radio background provides a key constraint on the 
counts of sources too faint to be individually detected. A 
summary of the available estimates of the background in- 
tensity as a function of frequency is presented in Fig. [5] 
where the dashed line shows the contribution of extragalac- 
tic sources computed using the present model at frequencies 
of up to 5 GHz and the De Zotti et al. (2005) model at 
higher frequencies. Our estimates agree very well with those 
by Gervasi et al. (2008), obtained fitting smooth functions 
to the observed source counts. 

The directly observed counts already account for most 
of the background intensity estimated by Bridle et al. (1967) 
and Wall et al. (1970), implying that no much room is left 
for new source populations and for truly diffuse emissions. 
The main contribution comes from AGNs. The contribution 
of star-forming galaxies increases with decreasing frequency, 
reaching a fraction of 39.5% percent at 150 MHz. 

On the other hand, after subtracting a model for the 
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AGN: Mauch & Sadler (2007) o 
SB: Mouch & Sadler 2007) X 
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Figure 7. Local luminosity functions at 1.4 GHz of radio AGNs 
(diamonds) and of star- forming galaxies (asterisks), as estimated 
by Mauch & Sadler (2007). The lines show the sum of the ex- 
trapolations to 1.4 GHz of the 60^tm local luminosity functions 
of "warm" (usually interpreted as starburst) and "cold" (normal 
late type) IRAS galaxies by Takeuchi et al. (2003): the dashed 
line refers to the linear radio/far-IR relationship of eq. JHJ, while 
the solid line is based on that of eq. J9j, which deviates from 
linearity at low luminosities. 
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Figure 8. Contributions of different populations of star-forming 
galaxies and of AGNs to the (sub)-mjy source counts at 1.4 GHz. 



Galactic emission and the CMB from the measurements of 
the second-generation balloon-borne experiment ARCADE- 
2, Fixsen et al. (2009) found excess radiation at 3 GHz about 
5 times brighter than the estimated contribution from ex- 
tragalactic radio sources. From a re-analysis of several large- 
area surveys at lower frequencies to separate the Galactic 
and extragalactic components, they derived background in- 
tensities much higher than yielded by earlier estimates. It 
may be noted, however, that their power-law fit to the back- 
ground spectrum exceeds the minimum total sky brightness 
temperature (including Galactic emission) measured by Bri- 
dle et al. (1967) at 81.5 MHz. 



Figure 9. Contributions of AGNs and star-forming galaxies to 
the extragalactic radio background as a function of frequency as 
yielded by the present model compared with observational esti- 
mates. Triangles and crosses refer to estimates based on the T-T 
plot method (Turtle et al. 1962), while squares refer to estimates 
by Fixsen et al. (2009) either using ARCADE 2 data (squares) or 
re-analyzing published data from large area surveys ( X ) . 



5 IMPLICATIONS FOR THE 

INTERPRETATION OF DATA ON 
LARGE-SCALE STRUCTURE 

Extragalactic radio sources are well suited to probe the 
large-scale structure of the Universe: detectable up to high 
redshifts, they are unaffected by dust extinction, and can 
thus provide an unbiased sampling of volumes larger than 
those probed by optical surveys. 



5.1 The angular correlation function and its 
implications 

Thanks to the deep radio surveys covering large areas of the 
sky, such as the NVSS (Condon et al. 1998), the SUMSS- 
MGPS (Mauch et al. 2003, Murphy et al. 2007), the FIRST 
(Becker et al. 1995), the WENSS (Rengelink et al. 1997), 
the angular correlation function, w(9), of radio sources has 
been measured with high statistical significance up to large 
angular scales (Cress et al. 1996; Magliocchetti et al. 1998; 
Rengelink 1999; Blake & Wall 2002; Overzier et al. 2003; 
Blake et al. 2004), well beyond those reached by optical sur- 
veys (Padmanabhan et al. 2009). Amongst the cited surveys, 
the NVSS is characterized by the most extensive sky cover- 
age and can thus provide the best clustering statistics. The 
w(8) measured for NVSS sources brighter than 10 mjy, is 
well described by a power-law of slope —0.8 extending from 
~ 0.1 degrees up to scales of almost 10 degrees (Blake & 
Wall 2002). 

The angular two-point correlation function of a given 
class of objects is related to the spatial correlation function, 
£(r,z), by Limber's (1953) equation: 



(9) = / dzAf 2 (z) / d(8z)£[r(Sz,8),z 



dzM{z) 



(10) 



where AT(z) is the redshift distribution and r(Sz, 8) is the 
comoving spatial distance between two objects located at 
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Figure 10. Two-point angular correlation function of NVSS 
sources with 51.40112 ^ lOmJy as measured by Blake & Wall 
(2002, circles) compared with the model by Negrello et al. (2006, 
red thin curves) and with the present model (blue thick curves) 
fitting the redshift distribution by Brookes et al. (2008). The 
dashed curves include the contribution of a constant offset e = 
0.0001 to w(9) in order to account for the effect of possible spu- 
rious density gradients in the survey. 



redshifts z and z + 5z and separated by an angle 8 on the 
sky. On sufficiently large scales, where the clustering signal is 
produced by galaxies residing in distinct dark-matter halos, 
and under the assumption of a one-to-one correspondence 
between sources and their host halos, £(r, z) can be written 
as the product of the correlation function of dark matter, 
£dm, times the square of the bias parameter, 6 (Matarrese 
et al. 1997; Moscardini et al. 1998): 

Z(r,z) = b 2 (M cS ,z)&M(r,z). (11) 

Here M e s is the effective mass of the dark matter haloes in 
which the sources reside and b is derived according to the 
prescriptions of Sheth & Tormen (1999). 

Earlier analyses generally adopted the N(z) from Dun- 
lop & Peacock (1990) models, that were found to be incon- 
sistent with the observational determination by Brookes et 
al. (2008). We have therefore repeated the analysis with the 
J\T(z) yielded by the present model that provides a good 
smooth fit of the latter data. As pointed out by Negrello et 
al. (2006), in the standard hierarchical clustering scenario, 
for z £ 1 angular scales 8 Si 2° correspond to linear scales 
where the correlation function is negative, so that the ob- 
servation of a positive correlation function on large angu- 
lar scales entails strong constraints on the bias parameter 
or, equivalently, on the effective halo mass, at high red- 
shifts. This is because the amplitude of w(8) on smaller 
scales, where the contributions from all the relevant red- 
shifts are positive, requires the local value of the effective 
mass, M e g(z = 0), to be very large (Negrello et al. 2006 
obtained M eB (z = 0) ~ 10 15 M Q /h). For such values of M e g 
the bias factor increases steeply with z, blowing up the neg- 
ative contributions to w(8). If M c ff is that large and almost 
independent of z, as it is found to be in the case of optically- 



selected QSOs (Porciani et al. 2004; Croom et al. 2004), w(8) 
must be negative on scales above 2°-3°, contrary to obser- 
vations. 

Although the new redshift distribution peaks at a some- 
what lower z than predicted by the pure luminosity evolu- 
tion Dunlop & Peacock (1990) model adopted by Negrello et 
al. (2006), and therefore the negative contributions to w(8) 
are lower, we still find that M e g cannot keep constant up 
to z > 1. The observed w(8) can be accounted for if M e s 
decreases with increasing z mirroring the evolution with cos- 
mic time of the characteristic halo mass entering the non- 
linear regime. The best fit is obtained with a somewhat lower 
value for the local effective mass, M eS (z = 0) ~ lO 14 ' 5 M /h 
(Fig. OH). 

5.2 Integrated Sachs-Wolfe (ISW) effect 

The ISW effect describes the influence of the evolution of 
the gravitational potential in time-variable, linear, metric 
perturbations on CMB photons that traverse them. A non- 
zero effect arises if the gravitational potential decays, as in 
the case of an open universe when the effect of the space 
curvature is important, or when the dynamics of the uni- 
verse is dominated by dark energy. As first pointed out by 
Crittenden & Turok (1996) a promising way of probing the 
ISW effect is through correlating fluctuations in the Cosmic 
Microwave Background (CMB) with large-scale structure. 
Since the timescale for the decay of the potential is of the 
order of the present-day Hubble time, the effect is largely 
canceled on small scales, because photons travel through 
multiple density peaks and troughs. This is why surveys 
covering large areas of the sky and probing the large scale 
distribution up to 2 ~ 1 are necessary. These are the prop- 
erties of the NVSS survey, which is therefore a particularly 
well-suited ISW probe, and indeed has been extensively ex- 
ploited to look for the ISW signal (Boughn & Crittenden 
2004, 2005; Pietrobon et al. 2006; McEwen et al. 2007, 2008; 
Ho et al. 2008; Giannantonio et al. 2008; Raccanelli et al. 
2008). 

The comparison of the correlations inferred from the 
data with model predictions requires once again the redshift 
distribution and the bias parameter as a function of red- 
shift. Again, all analyses carried out so far have used red- 
shift distributions inconsistent with the CENSORS results, 
and this has motivated our re-analysis exploiting the red- 
shift distribution yielded by the model. In addition we have 
redone the cross-correlation between the spatial distribution 
of NVSS sources brighter than 10 mjy and the CMB map 
from WMAP 5yr data provided by Delabrouille et al. (2009). 
We have also exploited the WMAP 5yr ILC map (Hinshaw 
et al. 2009) for cross-checks only, because it is declared to 
be sub-optimal for the angular scales we are interested in. 

Since foreground cleaning of CMB maps is never per- 
fect, we need to use a mask to remove the most contaminated 
areas. The choice of the mask must be done carefully. In fact, 
while a small mask allows a wider statistics to compute the 
CCF, it can also lead to a higher residual foreground con- 
tamination on the CMB map, thus diluting the CCF signal. 
We tried both the KQ85 and KQ75 masks provided by the 
WMAP team, which mask 18% and 28% of the sky respec- 
tively, including the strongest point sources in the WMAP 
channels. 
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The estimate of the NVSS/CMB correlation, as a func- 
tion of the angular scale, and of the associated errors, was 
performed as in Raccanelli et al. (2008), where the formal- 
ism used to compute the theoretical cross-correlation func- 
tion is also presented. The results for both CMB maps and 
masks are shown in the upper panel of Fig. [11] The CCFs 
computed using the extended KQ75 mask are significantly 
higher, suggesting that the KQ85 is not sufficient to avoid 
a significant effect of the residual contamination present in 
the CMB maps. Once the same mask is adopted, the curves 
corresponding to the different CMB maps are quite similar. 

In the following, we used the Delabrouille et al. (2009) 
CMB map with the extended KQ75 mask. To estimate the 
errors on the empirical CCF we simulated ~ 3800 NVSS 
maps by randomly redistributing the unmasked pixels of 
the true NVSS map. For each simulated NVSS map we com- 
puted the CCF with the CMB map. In this way we obtained 
a model-independent test of the null-hypothesis (i.e. no ISW 
signal). Since the distribution of sources in the simulated 
maps does not trace the large scale density field with which 
the CMB interacts, the distribution of the simulated val- 
ues of the CMB-NVSS cross-correlation function (CCF(6»)) 
allows us to estimate the level of significance at which the 
null hypothesis can be rejected. The results are shown in the 
lower panel of Fig. [11] The CCF differs from zero at 2.8 a 
for the first point and at 2.4-1.7 a for the following 4 points. 
Since the data points and their errors are correlated, the 
overall significance of the ISW signal cannot be straightfor- 
wardly derived from the significance of each point. There- 
fore, we have estimated it from the fraction of simulated 
maps yielding a CCF(#) at the observed level, or higher. 
This fraction turned out to be ~ 2.6 ■ 10~ 4 , which implies a 
significance of the ISW signal one order of magnitude higher 
than found in Raccanelli et al. (2008) and in other previous 
analyses. This improvement is due to a better CMB map, 
as allowed by 5 years of integration by WMAP, and, even 
more, to the use of the wider KQ75 mask. 

As shown by Fig. 111! lower panel, the cross-correlation 
power spectrum between the surface density fluctuations of 
NVSS sources and the CMB fluctuations expected for the 
"concordance" ACDM cosmology turns out to be in good 
agreement with the empirical determination. As in the case 
of w(ff) since we have now a reliable observational estimate 
of the redshift distribution of NVSS sources, the critical in- 
gredient is the evolution with redshift of the bias factor b(z) 
of radio sources. Without a reliable theory for it any strong 
inference from ISW estimates on the need of new gravita- 
tional physics on cosmological scales (Afshordi et al. 2008) 
or of a large local primordial non-Gaussianity (Afshordi & 
Tolley 2008) is premature. 



6 CONCLUSIONS 

We have worked out a new analytical evolutionary model 
that describes the population properties of radio sources at 
frequencies ^ 5 GHz, thus complementing the De Zotti et 

al. (2005) model, holding at higher frequencies. The model 
provides a satisfactory fit to a broad variety of data includ- 
ing luminosity functions of steep- and flat-spectrum sources, 
source counts at all frequencies from 0.15 to 5 GHz and over 
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Figure 11. WMAP-NVSS cross-correlation function (CCF). The 
upper panel compares CCF estimates obtained using different re- 
constructions of CMB maps from WMAP 5yr data (Delabrouille 
et al. 2009, dotted line; Hinshaw et al. 2009, dashed line) and dif- 
ferent masks (KQ85, diamonds, and KQ75, asterisks; masks pro- 
vided by the WMAP team). The solid line shows the expected 
signal for a concordance cosmology, based on the redshift distri- 
bution yielded by the present model and the bias function de- 
scribed in the text. In the lower panel, such signal is compared 
with our best estimate of the WMAP-NVSS CCF (CMB map by 
Delabrouille et al. 2009 and large mask; points with with la error 
bars) 



the full flux density range covered by observations, and red- 
shift distributions. 

We find that luminosity evolution is still sufficient to fit 
the wealth of available data on radio AGNs. However, the fit 
requires a luminosity-dependent decline of space densities of 
steep-spectrum sources at high redshifts, with positive evo- 
lution continuing up to higher redshifts for the more lumi- 
nous sources, thus confirming earlier indications of a "down- 
sizing" also for radio sources. The data do not require a simi- 
lar behaviour for FSRQs. However, the redshift at which the 
evolution peaks is found to be substantially higher for these 
sources than for the lower luminosity BL Lac objects, again 
consistent with "downsizing". Note that, while earlier stud- 
ies (Hook et al. 1998; Waddington et al. 2001; Vigotti et 
al. 2003; Cirasuolo et al. 2005, 2006) were based on limited 
data sets, the evidence presented here comes from a global 
analysis of all the available data. 

The upturn of source counts at sub-mjy levels is ac- 
counted for by a straightforward extrapolation, exploiting 
the well established correlation between far-IR and radio 
luminosities, of evolutionary models accounting for the far- 
IR counts and redshift distributions of star-forming galax- 
ies. As shown by earlier analyses (Yun et al. 2001; Best et 
al. 2005), a comparison of radio and 60^im local luminosity 
functions indicates deviations from a linear relation at low 
luminosities. We present a simple analytic expression de- 
scribing the correlation over the full luminosity range. The 
currently available surveys are not deep enough to provide 
an independent test of the deviations from linearity, but sur- 
veys with the new generation of radio telescopes (e.g. the 
SKA and its pathfinders) will. 

The directly counted sources already account for most 
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of the extragalactic background intensity determined by 
means of the classical T-T plot method. The dominant con- 
tribution comes from radio AGNs; the contribution of star- 
forming galaxies is higher at the lowest frequencies where it 
reaches 39.5% percent. 

We also discuss the implications of the new model for 
the interpretation of data on large-scale clustering of radio 
sources and on the Integrated Sachs-Wolfe effect, based pri- 
marily on the NVSS survey. Since both the sky coverage 
and the depth of this survey are much higher than those of 
large surveys in all the other wavebands, it provides unique 
constraints on the high-z evolution of the bias factor, b(z), 
on very large scales. The present analysis confirms and up- 
dates the conclusion of Negrello et al. (2006) that, if wc 
adopt the Sheth & Tormen (1999) formalism to derive b(z) 
as a function of the effective halo mass, M e g, associated 
to radio sources, M e g must decrease at z > 1 rather than 
keeping essentially constant as found for optically selected 
QSOs. In the absence of a reliable theory for the evolution of 
the bias parameter of radio sources, any claim on the need 
of new gravitational physics on cosmological scales or of a 
large local primordial non-Gaussianity, based on the spatial 
distribution of NVSS sources is, at best, premature. 
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